Purpose: To measure the capacity for each strain to grow as a biofilm

Protocol:

library(tidyverse)
library(xlsx)
library(broom)
library(rstatix)
library(ape)
library(ggbeeswarm)
library(formattable)
library(ggpubr)
library(cowplot)
library(corrplot)
library(ggpubr)
library(grid)
library(gridGraphics)

`%!in%` <- negate(`%in%`)

Formatting and set up

Short cuts:

strains <- c("14018", "315-A", "C0011E4", "JCP7275", "C0084H9", "JCP8108", "JCP8017B", "C0040C2" , "UM35",  "JCP8066", "C0093B3", "AMD", "C0096A1", "C0102A2",  "UM224", "Gv5-1", "C0179B4", "C0056B5", "Gv101", "C0100B2", "C0101A1", "CMW7778B")

strains0 <- c("14018", "315-A", "C0011E4", "JCP7275", "C0084H9", "JCP8108", "JCP8017B", "C0040C2" , "UM35",  "JCP8066", "C0093B3", "AMD", "C0096A1", "C0102A2",  "UM224", "5-1", "C0179B4", "C0056B5", "101", "C0100B2", "C0101A1", "CMW7778B")

species <- c("Gardnerella vaginalis", "Gardnerella sp. 2", "Gardnerella sp. 3", "Gardnerella piotii", "Gardnerella leopoldii", "Gardnerella swidsinskii", "Gardnerella sp. 7", "Gardnerella sp. 8", "Gardnerella sp. 10", "Gardnerella sp. 11", "Gardnerella sp. 12")

figureOut <- "../../experiments_figures"
suppTableOut <- "../../experiments_figures/supplementary_posthoc_tables.xlsx"

Import data

midlog_raw_data <- read_csv("./data/biofilm_midlog_raw_data.csv")
rawPlanktonicODs <- read_csv("./data/planktonic_od_raw_data.csv")
rawBiofilmODs <- read_csv("./data/biofilm_od_raw_data.csv")
rawSafPlanktonicODs <- read_csv("./data/safranin_planktonic_raw_data.csv")
rawSafBiofilmODs <- read_csv("./data/safranin_biofilm_raw_data.csv")

Plot settings

theme_set(theme_bw()+
          theme(panel.grid = element_blank()))

speciesColor <- scale_color_manual(values=c("#CC79A7", "#D55E00", "#56B4E9", "#0072B2", "#E69F00", "#009E73", 
                                                 "darkgreen", "gold3", "darkred", "darkblue",
                                                 "mediumpurple4"))

speciesFill <- scale_fill_manual(values=c("#CC79A7", "#D55E00", "#56B4E9", "#0072B2", "#E69F00", "#009E73", 
                                                 "darkgreen", "gold3", "darkred", "darkblue",
                                                 "mediumpurple4"))

Set up growth

Propagate twice and then incubate to mid log. 1:10 Dilution for mid Log: Protocol 1: 4 hour incubation Protocol 2: 7 hour incubation Protocol 3: 15 hour incubation ## Clean data

setupDF <- midlog_raw_data %>%
  group_by(Number, BioRep, Species, Strain, Batch, Step) %>%
  mutate(OD=OD-mean(.$OD[.$Number=="Blank"]),
                      OD=case_when(OD>0~OD,
                      OD<=0~0)) %>%
  filter(Number!="Blank") %>%
  mutate(Step=factor(Step, levels=c("Overnight 1", "Overnight 2", "Mid Log"))) %>%
  mutate(Strain=factor(Strain, levels=strains),
         Species=factor(Species, levels=species),
         Sample=paste(Strain, Batch, BioRep, sep="."),
         BatchNum=str_extract(Batch, "(?<=[A|B])[0-9]"),
         BatchRep=paste(BatchNum, BioRep, sep="."),
         StrainN=as.numeric(Strain),
         rectFill=case_when(StrainN %% 2 == 0 ~ "A",
                           StrainN %% 2 != 0 ~ "B")) %>%
  ungroup %>%
  with_groups(c("Sample", "Number", "BioRep", "Batch", "BatchNum", "BatchRep", "StrainN", "rectFill", "Species", "Strain", "Protocol",  "Step"), summarize, ODm=mean(OD), ODsd=sd(OD))

Assess mid log growth for removing poor growers with OD600 < 0.05

setupDF %>%
  filter(Step=="Mid Log") %>%
  ggplot(aes(x=Strain, y=ODm)) +
  geom_point(aes(shape=BioRep)) +
  geom_hline(aes(yintercept=0.05), linetype=2, color="gray") +
  theme(axis.text.x = element_text(angle=45, hjust=1)) +
  labs(x=NULL, y=bquote(OD[600]), shape="Biological Replicate")

Create quality filter for OD600 > 0.05

midLogFail <- setupDF %>%
  filter(Step=="Mid Log", 
         ODm < 0.05) %>%
  .$Sample
  
midLogFail
##  [1] "14018.1A1.A"    "14018.1A1.B"    "14018.1A2.A"    "14018.1A2.B"   
##  [5] "AMD.2A1.B"      "C0096A1.2A2.A"  "C0096A1.2A3.A"  "C0100B2.1AB3.A"
##  [9] "C0100B2.1AB3.B" "C0102A2.1A1.B"  "C0102A2.1A2.B"  "CMW7778B.2A1.A"
## [13] "CMW7778B.2A1.B" "CMW7778B.2A2.A" "CMW7778B.2A2.B" "CMW7778B.2A3.A"
## [17] "Gv5-1.1AB3.A"   "Gv5-1.1AB3.B"   "Gv5-1.1B1.A"    "Gv5-1.1B1.B"   
## [21] "Gv5-1.1B2.B"    "JCP7275.1AB3.A" "JCP7275.1AB3.B" "JCP7275.1AB4.B"
## [25] "JCP7275.1B1.A"  "JCP7275.1B1.B"  "JCP7275.1B2.A"  "JCP7275.1B2.B" 
## [29] "JCP8066.1A2.A"  "JCP8066.1A2.B"  "JCP8066.1AB3.A" "JCP8066.1AB3.B"
## [33] "JCP8108.1B1.A"  "JCP8108.1B2.A"  "JCP8108.1B2.B"  "UM224.1A1.A"   
## [37] "UM224.1A2.B"    "UM35.3A1.B"     "UM35.3A2.A"     "UM35.3A2.B"    
## [41] "UM35.3A3.B"
overnightPlot <- setupDF %>%
  filter(Sample %!in% midLogFail) %>%
  ggplot(aes(x=Strain, y=ODm)) +
  geom_point(aes(shape=BatchNum), alpha=0) +
  geom_rect(aes(xmin=StrainN-.5, xmax=StrainN+.5, ymin = 0, ymax = 1, fill=rectFill), alpha = 0.2, show.legend = FALSE) +
  geom_pointrange(aes(ymin=(ODm-ODsd), ymax=(ODm+ODsd), shape=BatchNum, color=Species), position = position_quasirandom(), size=0.3, alpha=0.7, show.legend = FALSE) +
  speciesColor +
  scale_fill_manual(values=c("white", "gray"), labels=c("A","B")) +
  ylim(0,1) +
  facet_grid(.~Step) +
  coord_flip() +
  scale_x_discrete(limits=rev) +
  theme(axis.text.x = element_text(size=8),
        axis.text.y = element_text(size=8),
        axis.title = element_text(size=11),
        panel.grid = element_blank(),
        legend.position = "none") +
  labs(x=NULL, y=bquote(OD[600]), shape="Experiment Batch")

#annotation bars
protocolBar <- setupDF %>%
  select(Strain, Species, Protocol) %>%
  mutate(Protocol=as.character(Protocol)) %>%
  unique %>%
  ggplot() +
    geom_bar(aes(x=Strain, y=2, fill=Protocol), stat="identity", width = 1.1) +
    scale_x_discrete(limits=rev) +
    scale_fill_manual(values=c("cornflowerblue", "burlywood1", "forestgreen")) +
    coord_flip() +
    theme_void() +
    theme(legend.position = "none")

speciesBar <- setupDF %>%
  select(Strain, Species, Protocol) %>%
  unique %>%
  ggplot() +
    geom_bar(aes(x=Strain, y=1, fill=Species), stat="identity", width=1.1) +
    scale_x_discrete(limits=rev) +
    speciesFill +
    coord_flip() +
    theme_void() +
    theme(legend.position = "none")

#Legends
protocol_legend <- ggpubr::get_legend(protocolBar +
    guides(color=guide_legend(ncol = 1)) +
    theme(legend.position = "bottom",
          legend.direction = "horizontal"))

species_legend <- ggpubr::get_legend(speciesBar +
    guides(color=guide_legend(ncol = 1)) +
    theme(legend.position = "bottom",
          legend.direction = "horizontal"))

experiment_legend <- ggpubr::get_legend(
  overnightPlot +
  guides(shape=guide_legend(nrow = 1, override.aes = list(alpha=1, size=4))) +
  theme(legend.position = "bottom"))


a <- plot_grid(overnightPlot, speciesBar, protocolBar, nrow = 1, align = "h", axis = "bt", rel_widths = c(1, 0.02, 0.02))
b <- plot_grid(a, as_ggplot(species_legend), nrow=2, rel_heights = c(1, 0.3))
c <- plot_grid(as_ggplot(protocol_legend), as_ggplot(experiment_legend), nrow=1)
plot_grid(b, c, ncol=1, nrow=2, rel_heights = c(1, 0.1))

Percent Biofilm

Calculate percent OD

planktonicODs <- rawPlanktonicODs %>%
  group_by(Batch) %>%
  mutate(planktonicOD=OD-mean(.$OD[.$Number=="Blank"]),
         planktonicOD=case_when(planktonicOD>0~planktonicOD,
                                planktonicOD<=0~0)) %>%
  filter(Number!="Blank") %>%
  select(-OD, -Measure) %>%
  ungroup

biofilmODs <- rawBiofilmODs %>%
  group_by(Batch) %>%
  mutate(biofilmOD=OD-mean(.$OD[.$Number=="Blank"]),
         biofilmOD=case_when(biofilmOD>0~biofilmOD,
                             biofilmOD<=0~0)) %>%
  filter(Number!="Blank") %>%
  select(-OD, -Measure) %>%
  ungroup

percentOD <- planktonicODs %>%
  left_join(biofilmODs) %>%
  mutate(Sample=paste(Strain, Batch, BioRep, sep=".")) %>%
  filter(Sample %!in% midLogFail) %>% # quality filter for poor growth at mid log phase 
  mutate(edgeWell=str_detect(well, pattern="[A|H]"),
         edgeWell=factor(edgeWell, levels=c(TRUE, FALSE)),
         innerWell=str_detect(well, pattern="[D|E]"),
         innerWell=factor(innerWell, levels=c(TRUE, FALSE)),
         topWell=str_detect(well, pattern="A"),
         topWell=factor(topWell, levels=c(TRUE, FALSE)))

head(percentOD)
## # A tibble: 6 × 13
##   well  Number BioRep Strain  Species      Protocol Batch planktonicOD biofilmOD
##   <chr> <chr>  <chr>  <chr>   <chr>           <dbl> <chr>        <dbl>     <dbl>
## 1 A2    1      A      C0011E4 Gardnerella…        1 1A1         0.196     0.0663
## 2 B2    1      A      C0011E4 Gardnerella…        1 1A1         0.193     0.0773
## 3 C2    1      A      C0011E4 Gardnerella…        1 1A1         0.172     0.0833
## 4 D2    1      A      C0011E4 Gardnerella…        1 1A1         0.138     0.0913
## 5 E2    1      B      C0011E4 Gardnerella…        1 1A1         0.0895    0.152 
## 6 F2    1      B      C0011E4 Gardnerella…        1 1A1         0.0845    0.221 
## # ℹ 4 more variables: Sample <chr>, edgeWell <fct>, innerWell <fct>,
## #   topWell <fct>
edgeWellPlot <- percentOD %>%
  mutate(percentBiofilm=(biofilmOD/(biofilmOD+planktonicOD))*100,
         Strain=factor(Strain, levels = strains),
         Species=factor(Species, levels = species)) %>%
  ggplot(aes(x=Strain, y=percentBiofilm, shape=BioRep, color=edgeWell)) +
  geom_quasirandom() + 
  theme(axis.text.x = element_text(angle=45, hjust=1), 
        legend.position = "bottom",
        legend.direction = "vertical") +
  labs(x=NULL, y="Percent Biofilm", shape="Biological Replicate", title="Edge Wells")

innerWellPlot <- percentOD %>%
  mutate(percentBiofilm=(biofilmOD/(biofilmOD+planktonicOD))*100,
         Strain=factor(Strain, levels = strains),
         Species=factor(Species, levels = species)) %>%
  ggplot(aes(x=Strain, y=percentBiofilm, shape=BioRep, color=innerWell)) +
  geom_quasirandom() + 
  theme(axis.text.x = element_text(angle=45, hjust=1),
        legend.position = "bottom",
        legend.direction = "vertical") +
  labs(x=NULL, y=NULL, shape="Biological Replicate", title="Inner Wells")

plot_grid(edgeWellPlot, innerWellPlot, nrow=1)

percentODDF0 <- percentOD %>%
  mutate(percentBiofilm=(biofilmOD/(biofilmOD+planktonicOD))*100,
         Strain=factor(Strain, levels = strains),
         Species=factor(Species, levels = species), 
         Protocol=as.character(Protocol), 
         BatchRep=paste(Batch, BioRep, sep=".")) %>%
  gather("measure", "value", c(planktonicOD, biofilmOD, percentBiofilm)) %>%
  mutate(measure=factor(measure, levels=c("biofilmOD", "planktonicOD", "percentBiofilm"), labels=c("Biofilm OD", "Planktonic OD", "Percent OD")),
         BatchNum=str_extract(Batch, "(?<=[A|B])[0-9]"),
         BatchRep=paste(BatchNum, BioRep, sep="."))

percentODDF0 %>%
  ggplot(aes(x=Strain, y=value)) +
    geom_boxplot(alpha=0) +
    geom_quasirandom(aes(color=edgeWell), alpha=0.7) +
    scale_y_continuous(position = "right") +
    facet_grid(measure~., scales="free") +
    theme(axis.text.x = element_text(angle=45, hjust=1)) +
    labs(x=NULL, y=NULL, shape="Biological Replicate")

Remove Edge Wells

naSamples <- percentODDF0 %>%
  filter(edgeWell==FALSE) %>%
  spread(measure, value) %>%
  filter(is.na(`Percent OD`)) %>%
  .$Sample %>%
  unique

naSamples
## [1] "C0093B3.3A2.A" "C0093B3.3A2.B"
percentODDF <- percentODDF0 %>%
  filter(edgeWell==FALSE, # remove samples
         Sample %!in% naSamples) %>% # remove samples with no growth
  mutate(method="Percent OD")

Isolates with extra samples

non4StrainsPOD <- percentODDF0 %>%
 mutate(Sample=paste(Number, Sample, sep=".")) %>%
 with_groups(Strain, summarize, n=n_distinct(Sample)) %>%
 filter(n!=4) %>%
 .$Strain

percentODDF %>%
  filter(Strain %in% non4StrainsPOD) %>%
  mutate(Sample=paste(Number, Sample, sep=".")) %>%
  with_groups(Strain, mutate, n=n_distinct(Sample)) %>%
  with_groups(c("Number", "Strain", "Species", "Protocol", "Batch", "method", "Sample", "BatchRep", "measure", "BatchNum", "n"), summarize, valueM=mean(value), valueSD=sd(value)) %>%
  mutate(labelYpos=case_when(measure=="Biofilm OD"~0.75,
                             measure=="Planktonic OD"~0.75,
                             measure=="Percent OD"~110)) %>%
  ggplot() +
  geom_boxplot(aes(x=Strain, y=valueM), alpha=0) +
  geom_pointrange(aes(x=Strain, y=valueM, ymin=valueM-valueSD, ymax=valueM+valueSD, shape=BatchRep), alpha=0.7, size=0.25, position=position_quasirandom()) + 
  geom_text(aes(x=Strain, y=labelYpos, label=n)) +
  scale_y_continuous(position = "right") +
  facet_grid(measure~., scales="free") +
  theme(axis.text.x = element_text(angle=45, hjust=1)) +
  labs(x=NULL, y=NULL)

Plot

## annotation portions of plot
#protocol bar
protocolBar <- percentODDF %>%
  select(Strain, Species, Protocol) %>%
  unique %>%
  ggplot() +
    geom_bar(aes(x=Strain, y=2, fill=Protocol), stat="identity", width = 1.1) +
    scale_fill_manual(values=c("cornflowerblue", "burlywood1", "forestgreen")) +
    theme_void() +
    theme(legend.position = "none")

# species bar
speciesBar <- percentODDF %>%
  select(Strain, Species, Protocol) %>%
  unique %>%
  ggplot() +
    geom_bar(aes(x=Strain, y=1, fill=Species), stat="identity", width=1.1) +
    speciesFill +
    theme_void() +
    theme(legend.position = "none")

percentODPlot <- percentODDF %>%
  with_groups(c("Number", "Strain", "Species", "Protocol", "Batch", "method", "Sample", "BatchRep", "measure", "BatchNum"), summarize, valueM=mean(value), valueSD=sd(value)) %>%
  ggplot() +
  geom_boxplot(aes(x=Strain, y=valueM), alpha=0) +
  geom_pointrange(aes(x=Strain, y=valueM, ymin=valueM-valueSD, ymax=valueM+valueSD, color=Species), alpha=0.7, size=0.25, position=position_quasirandom()) + 
  speciesColor +
  facet_grid(measure~., scales="free", switch = "y") +
  theme(axis.text.x = element_text(angle=45, hjust=1),
        legend.position = "none",
        strip.placement = "outside") +
  labs(x=NULL, y=NULL)

ODPlots <- plot_grid(protocolBar, speciesBar, percentODPlot, ncol = 1, rel_heights = c(0.06, 0.06, 1), align = "v", axis = "lr")
ODPlots

Compare Percent OD

Isolate ANOVA

strain_OD_DF <- percentODDF %>%
  spread(measure, value) %>%
  select(-`Biofilm OD`, -`Planktonic OD`) %>%
  dplyr::rename(percentOD=`Percent OD`) %>%
  with_groups(c(Sample, Species, Strain, Batch, BioRep), summarize, mPercentOD=mean(percentOD), sdPercentOD=sd(percentOD))

strain_OD_DF %>%
  ggplot(aes(sample=mPercentOD)) +
  geom_qq_line() +
  geom_qq()

#### One-way ANOVA

strain_OD_DF %>%
  anova_test(mPercentOD~Strain) %>%
  formattable
Effect DFn DFd F p p<.05 ges
Strain 20 55 5.49 2.6e-07
0.666

Tukey’s HSD

ODPostHoc <- strain_OD_DF %>%
  tukey_hsd(mPercentOD~Strain, p.adjust.method = "BH")

# ODPostHoc %>%
#   arrange(p.adj) %>%
#   write.xlsx(file=suppTableOut, sheetName="TableS5_percent_od_posthoc", col.names=TRUE, append=TRUE)
pvals <- ODPostHoc %>%
  select(group1, group2, p.adj) %>%
  dplyr::rename(var1=group2, var2=group1, cor=p.adj) %>%
  cor_spread() %>%
  column_to_rownames(var="rowname") %>%
  as.matrix()

ODPostHoc %>%
  select(group1, group2, estimate) %>%
  dplyr::rename(var1=group2, var2=group1, cor=estimate) %>%
  cor_spread() %>%
  column_to_rownames(var="rowname") %>%
  as.matrix() %>%
  corrplot(is.corr=FALSE, type = "lower", tl.col="black", p.mat = pvals, insig = "label_sig", sig.level = c(0.001, 0.01, 0.05), pch.cex = 0.9)

grid.echo()

percODPostHocPlot <- grid.grab()

Percent biofilm of Strains

strain_OD_DF %>%
  with_groups(Strain, summarize, `Percent OD Biofilm`=mean(mPercentOD)) %>%
  arrange(-`Percent OD Biofilm`) %>%
  formattable()
Strain Percent OD Biofilm
C0040C2 97.58907
JCP8066 91.16070
C0056B5 90.17673
C0096A1 84.37130
C0100B2 84.26520
C0102A2 83.82580
C0179B4 81.63978
Gv5-1 80.78987
UM224 80.75329
AMD 76.27148
UM35 70.98144
C0101A1 69.57704
JCP8017B 64.14116
315-A 61.51312
JCP8108 59.02140
C0011E4 51.01499
C0084H9 50.11188
Gv101 48.29067
C0093B3 41.46202
JCP7275 36.34031
CMW7778B 27.52904

Average of all strains

strain_OD_DF %>%
  summarize(`Average Percent OD Biofilm`=mean(mPercentOD)) %>%
  formattable()
Average Percent OD Biofilm
70.12149

Moran’s I

Load tree

We choose the weights as \(_{wij} = 1/d_{ij}\) , where the d’s is the distances measured on the tree:

We can now perform the analysis with Moran’s I:

Safranin

safPlankODDF <- rawSafPlanktonicODs %>% 
  group_by(Batch) %>%
  mutate(safPlanktonicOD=OD-mean(.$OD[.$Number=="Blank"]),
                      OD=case_when(OD>0~OD,
                               OD<=0~0)) %>%
  filter(Number!="Blank") %>%
  ungroup %>%
  mutate(row=str_extract(well, pattern="[A-H]")) %>%
  select(-OD, -Measure, -well)

safSafODDF <- rawSafBiofilmODs %>% 
  group_by(Batch) %>%
  mutate(safraninOD=OD-mean(.$OD[.$Number=="Blank"]),
                      OD=case_when(OD>0~OD,
                               OD<=0~0)) %>%
  filter(Number!="Blank") %>%
  ungroup %>%
  mutate(row=str_extract(well, pattern="[A-H]")) %>%
  select(-OD, -Measure, -well)


safraninOD1 <- safPlankODDF %>%
  full_join(safSafODDF) %>%
  select(row, everything(), safPlanktonicOD, safraninOD) %>%
  mutate(Sample=paste(Strain, Batch, BioRep, sep=".")) %>%
  filter(Sample %!in% midLogFail) %>% # quality filter for poor growth at mid log phase 
  mutate(edgeWell=str_detect(row, pattern="[A|H]"),
         edgeWell=factor(edgeWell, levels=c(TRUE, FALSE)),
         innerWell=str_detect(row, pattern="[D|E]"),
         innerWell=factor(innerWell, levels=c(TRUE, FALSE)),
         topWell=str_detect(row, pattern="A"),
         topWell=factor(topWell, levels=c(TRUE, FALSE)))
head(safraninOD1)
## # A tibble: 6 × 13
##   row   Number BioRep Strain  Species  Protocol Batch safPlanktonicOD safraninOD
##   <chr> <chr>  <chr>  <chr>   <chr>       <dbl> <chr>           <dbl>      <dbl>
## 1 A     1      A      C0011E4 Gardner…        1 1A1             0.151      0.202
## 2 B     1      A      C0011E4 Gardner…        1 1A1             0.142      0.217
## 3 C     1      A      C0011E4 Gardner…        1 1A1             0.136      0.232
## 4 D     1      A      C0011E4 Gardner…        1 1A1             0.159      0.201
## 5 E     1      B      C0011E4 Gardner…        1 1A1             0.152      0.179
## 6 F     1      B      C0011E4 Gardner…        1 1A1             0.116      0.103
## # ℹ 4 more variables: Sample <chr>, edgeWell <fct>, innerWell <fct>,
## #   topWell <fct>

Organize data and create dataframe

safDF <- safraninOD1 %>%
  mutate(Strain=factor(Strain, levels = strains),
         Species=factor(Species, levels = species),
         Protocol=as.character(Protocol),
         BatchNum=str_extract(Batch, "(?<=[A|B])[0-9]"),
         BatchRep=paste(BatchNum, BioRep, sep="."), 
         method="Safranin Stain")  %>%
  filter(Sample %!in% naSamples) # remove samples with no growth

Check if edge wells should be removed

safDF %>%
  ggplot(aes(x=Strain, y=safraninOD)) +
    geom_boxplot(alpha=0) +
    geom_quasirandom(aes(color=edgeWell), alpha=0.7) +
    scale_y_continuous(position = "right") +
    theme(axis.text.x = element_text(angle=45, hjust=1)) +
    labs(x=NULL, y=NULL, shape="Biological Replicate")

Do not remove edge wells ## Extra biological replicates

non4StrainsSaf <- safDF %>%
 mutate(Sample=paste(Number, Sample, sep=".")) %>%
 with_groups(Strain, summarize, n=n_distinct(Sample)) %>%
 filter(n!=4) %>%
 .$Strain

safDF %>%
  filter(Strain %in% non4StrainsSaf) %>%
  mutate(Sample=paste(Number, Sample, sep=".")) %>%
  with_groups(Strain, mutate, n=n_distinct(Sample)) %>%
  with_groups(c("Number", "Strain", "Species", "Protocol", "Batch", "method", "Sample", "BatchRep", "BatchNum", "n"), summarize, safraninM=mean(safraninOD), safraninSD=sd(safraninOD)) %>%
  ggplot() +
  geom_boxplot(aes(x=Strain, y=safraninM), alpha=0) +
  geom_pointrange(aes(x=Strain, y=safraninM, ymin=safraninM-safraninSD, ymax=safraninM+safraninSD, shape=BatchRep), alpha=0.7, position=position_quasirandom(), size=0.3) + 
  geom_text(aes(x=Strain, y=0.575, label=n)) +
  theme(axis.text.x = element_text(angle=45, hjust=1)) +
  labs(x=NULL, y=bquote("Safranin OD"[450]), shape="Batch")

Plot

safPlot <- safDF %>%
  ggplot(aes(x=Strain, y=safraninOD, color=Species, shape=BatchRep)) +
  geom_quasirandom(alpha=0.75) + 
  theme(axis.text.x = element_text(angle=45, hjust=1)) +
  labs(x=NULL, y=bquote("Safranin OD"[450]), shape="Biological Replicate")

safPlot2 <- safDF %>%
  with_groups(c("Number", "Strain", "Species", "Protocol", "Batch", "method", "Sample", "BatchRep", "BatchNum"), summarize, safraninM=mean(safraninOD), safraninSD=sd(safraninOD)) %>%
  ggplot() +
  geom_boxplot(aes(x=Strain, y=safraninM), alpha=0) +
  geom_pointrange(aes(x=Strain, y=safraninM, ymin=safraninM-safraninSD, ymax=safraninM+safraninSD, color=Species), alpha=0.7, position=position_quasirandom(), size=0.3) + 
  speciesColor +
  theme(axis.text.x = element_text(angle=45, hjust=1),
        legend.position = "none") +
  labs(x=NULL, y=bquote("Safranin OD"[450]), shape="Batch")

safPlots <- plot_grid(protocolBar, speciesBar, safPlot2, ncol = 1, rel_heights = c(0.06, 0.06, 1), align = "v", axis = "lr")

safPlots

Compare Safranin

Isolate ANOVA

safDF %>%
  with_groups(c(Sample, Species, Strain, Batch, BioRep), summarize, mSafOD=mean(safraninOD), sdSafOD=sd(safraninOD)) %>%
  ggplot(aes(sample=mSafOD)) +
  geom_qq_line() +
  geom_qq()

#### One-way ANOVA

strain_saf_DF <- safDF %>%
  with_groups(c(Sample, Species, Strain, Batch, BioRep), summarize, mSafOD=mean(safraninOD), sdSafOD=sd(safraninOD))

strain_saf_DF %>%
  anova_test(mSafOD~Strain) %>%
  formattable
Effect DFn DFd F p p<.05 ges
Strain 20 55 4.469 5.31e-06
0.619

Tukey’s HSD

safPostHoc <- strain_saf_DF %>%
  tukey_hsd(mSafOD~Strain, p.adjust.method = "BH")

# safPostHoc %>%
#   arrange(p.adj) %>%
#   write.xlsx(file=suppTableOut, sheetName="TableS6_safranin_posthoc", col.names=TRUE, append=TRUE)
pvals <- safPostHoc %>%
  select(group1, group2, p.adj) %>%
  dplyr::rename(var1=group2, var2=group1, cor=p.adj) %>%
  cor_spread() %>%
  column_to_rownames(var="rowname") %>%
  as.matrix()

safPostHoc %>%
  select(group1, group2, estimate) %>%
  dplyr::rename(var1=group2, var2=group1, cor=estimate) %>%
  cor_spread() %>%
  column_to_rownames(var="rowname") %>%
  as.matrix() %>%
  corrplot(is.corr=FALSE, type = "lower", tl.col="black", p.mat = pvals, insig = "label_sig", sig.level = c(0.001, 0.01, 0.05), pch.cex = 0.9)

grid.echo()

safraninPostHocPlot <- grid.grab()

Moran’s I

We can now perform the analysis with Moran’s I:

Figure 4: Biofilm growth as percent OD and Safranin stain

plot_grid(protocolBar, speciesBar, percentODPlot + speciesColor + theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()), safPlot2, ncol = 1, rel_heights = c(0.07, 0.07, 1, 0.5), align = "v", axis = "lr", labels = c("A", "", "", "B"), label_size = 15)

#ggsave(file.path(figureOut, paste(Sys.Date(), "Figure4_strain_biofilm.png", sep="_")))

Figure S3: Post Hoc Test Results

plot_grid(ggplot()+geom_blank(), percODPostHocPlot, ggplot()+geom_blank(), safraninPostHocPlot, nrow = 1, rel_widths = c(0.1, 1, 0.1, 1), labels = c(" ", "A", " ", "B"), label_size = 20)

#ggsave(file.path(figureOut, paste(Sys.Date(), "FigureS3_biofilm_posthocs.png", sep = "_")))

Session Info

sessionInfo()
## R version 4.4.0 (2024-04-24)
## Platform: x86_64-apple-darwin20
## Running under: macOS Ventura 13.6.6
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: America/New_York
## tzcode source: internal
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] gridGraphics_0.5-1 corrplot_0.92      cowplot_1.1.3      ggpubr_0.6.0      
##  [5] formattable_0.2.1  ggbeeswarm_0.7.2   ape_5.8            rstatix_0.7.2     
##  [9] broom_1.0.6        xlsx_0.6.5         lubridate_1.9.3    forcats_1.0.0     
## [13] stringr_1.5.1      dplyr_1.1.4        purrr_1.0.2        readr_2.1.5       
## [17] tidyr_1.3.1        tibble_3.2.1       ggplot2_3.5.1      tidyverse_2.0.0   
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.5      beeswarm_0.4.0    xfun_0.46         bslib_0.7.0      
##  [5] htmlwidgets_1.6.4 rJava_1.0-11      lattice_0.22-6    tzdb_0.4.0       
##  [9] vctrs_0.6.5       tools_4.4.0       generics_0.1.3    parallel_4.4.0   
## [13] fansi_1.0.6       highr_0.11        pkgconfig_2.0.3   lifecycle_1.0.4  
## [17] farver_2.1.2      compiler_4.4.0    munsell_0.5.1     carData_3.0-5    
## [21] vipor_0.4.7       htmltools_0.5.8.1 sass_0.4.9        yaml_2.3.9       
## [25] crayon_1.5.3      pillar_1.9.0      car_3.1-2         jquerylib_0.1.4  
## [29] cachem_1.1.0      abind_1.4-5       nlme_3.1-165      tidyselect_1.2.1 
## [33] digest_0.6.36     stringi_1.8.4     labeling_0.4.3    fastmap_1.2.0    
## [37] colorspace_2.1-0  cli_3.6.3         magrittr_2.0.3    xlsxjars_0.6.1   
## [41] utf8_1.2.4        withr_3.0.0       scales_1.3.0      backports_1.5.0  
## [45] bit64_4.0.5       timechange_0.3.0  rmarkdown_2.27    bit_4.0.5        
## [49] ggsignif_0.6.4    hms_1.1.3         evaluate_0.24.0   knitr_1.48       
## [53] rlang_1.1.4       Rcpp_1.0.13       glue_1.7.0        vroom_1.6.5      
## [57] rstudioapi_0.16.0 jsonlite_1.8.8    R6_2.5.1